knitr::opts_chunk$set(echo = TRUE)
library(st514)
library(car)
## Loading required package: carData
rm(list=ls())
options(scipen=999)

Opgave 5.1

I denne opgave skal vi evaluere Hotelling’s \(T^2\) og teste hypotesen \(H_0: \boldsymbol{\mu}' = [7, 11]\) baseret på et datasæt.

Opgave 5.1.a

Evaluer \(T^2\) for test af \(H_0: \boldsymbol{\mu}' = [7, 11]\) ved brug af data \[\mathbf{X} = \begin{bmatrix} 2 & 12 \\ 8 & 9 \\ 6 & 9 \\ 8 & 10 \end{bmatrix}\]

# Definer matricen X
x <- matrix(c(2, 8, 6, 8, 12, 9, 9, 10), 4, 2)

# Bestem stikprøvestørrelsen n og antal variable p
n <- nrow(x)
p <- ncol(x)

# Beregn stikprøvegennemsnittet
x.mean <- colMeans(x)

# Beregn kovariansmatricen
S <- cov(x)

# Definer middelværdivektoren under nulhypotesen
mu0 <- c(7, 11)

# Beregn Mahalanobis-afstanden i kvadrat
D2 <- t(x.mean - mu0) %*% solve(S) %*% (x.mean - mu0)

# Beregn Hotelling's T^2
T2 <- n * D2

# Vis resultaterne
cat("Stikprøvestørrelse n:", n, "\n",
    "Antal variable p:", p, "\n",
    "Stikprøvegennemsnit:", x.mean, "\n",
    "Middelværdivektor under H0:", mu0, "\n",
    "Kvadreret Mahalanobis-afstand D^2:", D2, "\n",
    "Hotelling's T^2:", T2, "\n")
## Stikprøvestørrelse n: 4 
##  Antal variable p: 2 
##  Stikprøvegennemsnit: 6 10 
##  Middelværdivektor under H0: 7 11 
##  Kvadreret Mahalanobis-afstand D^2: 3.409091 
##  Hotelling's T^2: 13.63636

Opgave 5.1.b

Angiv fordelingen af \(T^2\) for situationen i (a).

Hotelling’s \(T^2\) fordeling følger \(T^2 \sim \frac{(n-1)p}{n-p} \cdot F_{p, n-p}\), hvor: - \(n = 4\) (stikprøvestørrelse) - \(p = 2\) (antal variable)

Dette giver: \[T^2 \sim \frac{(4-1) \cdot 2}{4-2} \cdot F_{2, 2} = 3 \cdot F_{2, 2}\]

Under nulhypotesen \(H_0: \boldsymbol{\mu}' = [7, 11]\) følger \(T^2\) denne fordeling, og vi kan bruge dette til at beregne kritiske værdier og p-værdier.

Opgave 5.1.c

Brug (a) og (b) til at teste \(H_0\) på \(\alpha = .05\) niveau. Hvilken konklusion drager du?

# Beregn den kritiske værdi for T^2
alpha <- 0.05
critical_value <- 3 * qf(alpha, p, n - p, lower.tail = FALSE)

# Sammenlign T^2 med den kritiske værdi
reject_H0 <- T2 > critical_value

# Alternativt kan vi beregne p-værdien
Fval <- T2 * (n - p) / ((n - 1) * p)
pval <- pf(Fval, p, n - p, lower.tail = FALSE)

# Vis resultaterne
cat("Kritisk værdi for T^2 ved α =", alpha, ":", critical_value, "\n",
    "Beregnet T^2:", T2, "\n",
    "Forkast H0:", reject_H0, "\n",
    "F-værdi:", Fval, "\n",
    "p-værdi:", pval, "\n")
## Kritisk værdi for T^2 ved α = 0.05 : 57 
##  Beregnet T^2: 13.63636 
##  Forkast H0: FALSE 
##  F-værdi: 4.545455 
##  p-værdi: 0.1803279

Konklusion: Da den beregnede \(T^2\) er mindre end den kritiske værdi ved \(\alpha = 0.05\), kan vi ikke forkaste nulhypotesen \(H_0: \boldsymbol{\mu}' = [7, 11]\). Dette understøttes af p-værdien på 0.18, som er større end signifikansniveauet på 0.05. Der er derfor ikke statistisk belæg for at påstå, at populationsmiddelværdien er forskellig fra vektoren [7, 11].

Opgave 5.2

Brug data i Eksempel 5.1 og verificer at \(T^2\) forbliver uændret hvis hver observation \(\mathbf{x}_j, j = 1, 2, 3\) erstattes med \(\mathbf{C}\mathbf{x}_j\), hvor \[\mathbf{C} = \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}\]

Bemærk at observationerne \[\mathbf{C}\mathbf{x}_j = \begin{bmatrix} x_{j1} - x_{j2} \\ x_{j1} + x_{j2} \end{bmatrix}\] giver data-matricen \[\begin{bmatrix} (6 - 9) & (10 - 6) & (8 - 3) \\ (6 + 9) & (10 + 6) & (8 + 3) \end{bmatrix}'\]

# Definer oprindelige data fra Eksempel 5.1
x <- matrix(c(6, 10, 8, 9, 6, 3), 3, 2)

# Definer transformationsmatricen C
C <- matrix(c(1, 1, -1, 1), 2, 2)

# Transformer data ved at multiplicere hver observation med C
cx <- t(C %*% t(x))

# Vis original matrix
print(x)
##      [,1] [,2]
## [1,]    6    9
## [2,]   10    6
## [3,]    8    3
# Vis transformeret matrix
print(cx)
##      [,1] [,2]
## [1,]   -3   15
## [2,]    4   16
## [3,]    5   11
# Beregn T^2 for de originale data
n <- nrow(x)
p <- ncol(x)
x.mean <- colMeans(x)
S <- cov(x)
mu0 <- c(9, 5)
D2_original <- t(x.mean - mu0) %*% solve(S) %*% (x.mean - mu0)
T2_original <- n * D2_original

# Beregn T^2 for de transformerede data
n_transformed <- nrow(cx)
p_transformed <- ncol(cx)
cx.mean <- colMeans(cx)
S_transformed <- cov(cx)
mu0_transformed <- C %*% mu0
D2_transformed <- t(cx.mean - mu0_transformed) %*% solve(S_transformed) %*% (cx.mean - mu0_transformed)
T2_transformed <- n_transformed * D2_transformed

# Vis resultaterne
cat("\nHotelling's T^2 for originale data:", T2_original, "\n",
    "Hotelling's T^2 for transformerede data:", T2_transformed, "\n",
    "Forskel:", abs(T2_original - T2_transformed), "\n")
## 
## Hotelling's T^2 for originale data: 0.7777778 
##  Hotelling's T^2 for transformerede data: 0.7777778 
##  Forskel: 0

Konklusion: Resultaterne viser, at Hotelling’s \(T^2\) forbliver uændret under den lineære transformation med matricen \(\mathbf{C}\). Dette demonstrerer en vigtig egenskab ved Hotelling’s \(T^2\) statistik - den er invariant under lineære transformationer. Dette skyldes, at \(T^2\) måler den standardiserede kvadrerede Mahalanobis-afstand mellem stikprøvegennemsnittet og den hypoteserede middelværdi, og denne afstand bevares under lineære transformationer.

Opgave 5.3

I denne opgave skal vi evaluere \(T^2\) ved hjælp af alternative formler og sammenligne med Wilks’ lambda.

Opgave 5.3.a

Brug formel (5-15) til at evaluere \(T^2\) for data i Opgave 5.1.

# Data fra Opgave 5.1
x <- matrix(c(2, 8, 6, 8, 12, 9, 9, 10), 4, 2)
n <- nrow(x)
p <- ncol(x)
mu0 <- c(7, 11)

# Beregn S (SSCP-matricen - Sum of Squares and Cross Products)
# Bemærk at vi skal bruge (n-1)*S som er SSCP-matricen
S_sscp <- cov(x) * (n - 1)

# Beregn S0 (SSCP-matricen for afvigelser fra mu0)
S0 <- matrix(0, p, p)
for (i in 1:n) {
  deviation <- x[i,] - mu0
  S0 <- S0 + (deviation %*% t(deviation))
}

# Beregn T^2 ved brug af formel (5-15)
# (5-15): N-1 * (|S_0| / (|S_sscp|-1))
T2_alt <- (n-1) * (det(S0) / det(S_sscp) - 1)

# Sammenlign med den tidligere beregnede T^2
# Det er resultatet fra opgave 5.1.c
T2_original <- n * t(colMeans(x) - mu0) %*% solve(cov(x)) %*% (colMeans(x) - mu0)

# SSCP-matrix
print(S_sscp)
##      [,1] [,2]
## [1,]   24  -10
## [2,]  -10    6
# SSCP-matrix S0 (afvigelser fra mu0)
print(S0)
##      [,1] [,2]
## [1,]   28   -6
## [2,]   -6   10
# Vis resultaterne
cat("\nT^2 beregnet med formel (5-15):", T2_alt, "\n",
    "T^2 beregnet med standardformlen:", T2_original, "\n",
    "Forskel:", abs(T2_alt - T2_original), "\n")
## 
## T^2 beregnet med formel (5-15): 13.63636 
##  T^2 beregnet med standardformlen: 13.63636 
##  Forskel: 0.00000000000001953993

Opgave 5.3.b

Brug data i Opgave 5.1 til at evaluere \(\Lambda\) i (5-13). Evaluer også Wilks’ lambda.

# Beregn Lambda (formel 5-13)
lambda <- (det(S_sscp) / det(S0))^(n/2)

# Beregn Wilks' lambda
wilks_lambda <- det(S_sscp) / det(S0)

# Vis resultaterne
cat("Lambda (fra formel 5-13):", lambda, "\n",
    "Wilks' lambda:", wilks_lambda, "\n")
## Lambda (fra formel 5-13): 0.03251814 
##  Wilks' lambda: 0.1803279

Konklusion for opgave 5.3: Vi har evalueret Hotelling’s \(T^2\) ved hjælp af en alternativ formel (5-15) og verificeret, at den giver samme resultat som den oprindelige beregning. Derudover har vi beregnet \(\Lambda\) (udtryk 5-13) og Wilks’ lambda. Wilks’ lambda er et alternativt teststatistik, der ligesom \(T^2\) anvendes til at teste hypoteser om middelværdivektorer. En mindre værdi af Wilks’ lambda (tættere på 0) indikerer større afvigelse fra nulhypotesen.

Opgave 5.4

Brug sved-dataene i Tabel 5.1. (Se Eksempel 5.2.) Det er data om 20 sunde kvinder, hvor der bliver målt svedmængde, sodium i sveden og potassium i sveden.

Opgave 5.4.a

Bestem akserne på det 90% konfidensellipsoide for \(\boldsymbol{\mu}\). Bestem længderne af disse akser.

# Indlæs sved-data
data("T5-1")
colnames(tbl) <- c("Sved", "Sodium", "Potassium")

# Lav scatter plots for alle par af variable
pairs(tbl, main = "Scatter plots af sved-data")

# Beregn egenvektorer og egenværdier af kovariansmatricen
r <- eigen(cov(tbl))

# Antal observationer og variable
n <- nrow(tbl)
p <- ncol(tbl)

# Beregn den kritiske F-værdi for 90% konfidensniveau
Fcv <- qf(0.9, p, n - p)

# Beregn længden af akserne
axis_lengths <- rep(0, length(r$values))
for (i in 1:p) {
  axis_lengths[i] <- 2 * sqrt(r$values[i] * p * (n - 1) * Fcv / (n * (n - p)))
}

# Aksevektorer (egenvektorer af kovariansmatricen)
print(r$vectors)
##             [,1]        [,2]        [,3]
## [1,] -0.05084144 -0.57370364  0.81748351
## [2,] -0.99828352  0.05302042 -0.02487655
## [3,]  0.02907156  0.81734508  0.57541452
# Vis resultater
cat("Egenværdier:\n", r$values, "\n\n",
    "Akselængder for 90% konfidensellipsoide:\n", axis_lengths, "\n")
## Egenværdier:
##  200.4625 4.531591 1.301392 
## 
##  Akselængder for 90% konfidensellipsoide:
##  18.10135 2.721571 1.458473

Fortolkning af egenvektorer

  1. Den første aksevektor (kolonne 1 i egenvektor-matricen) viser den retning i rummet, hvor data har størst varians. Den udgør hovedaksen i konfidensellipsoiden.
  2. Den anden aksevektor (kolonne 2) er den næststørste variansretning, vinkelret på den første.
  3. Den tredje aksevektor (kolonne 3) er den mindste variansretning, vinkelret på både første og anden aksevektor.

Akselængderne repræsenterer diametrene af konfidensellipsoiden langs de respektive akser. Den største akselængde svarer til den retning, hvor der er størst usikkerhed om middelværdivektoren.

Opgave 5.4.b

Konstruer Q-Q plots for observationerne af svedraten, natrium-indhold og kaliumindhold. Konstruer de tre mulige scatter plots for par af observationer. Synes den multivariate normalfordelingsantagelse at være berettiget i dette tilfælde? Kommenter.

# Lav Q-Q plots for hver af de tre variable
par(mfrow = c(2, 2))

for (i in 1:p) {
  qqnorm(tbl[,i], main = colnames(tbl)[i])
  qqline(tbl[,i])
}

# Beregn Mahalanobis-afstandene
Sx <- cov(tbl)
D2 <- mahalanobis(tbl, colMeans(tbl), Sx)

# Lav et chi-i-anden Q-Q plot for multivariate normalitet
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

# Reset til 1x1 plotting (for fremtidige plots)
par(mfrow = c(1, 1))

Konklusion: Baseret på Q-Q plots for de enkelte variable ser vi, at dataene følger normalfordelingen ret godt, med kun mindre afvigelser i halerne. chi-i-anden Q-Q plottet for Mahalanobis-afstandene viser også en rimelig god tilpasning til en ret linje, hvilket understøtter antagelsen om multivariat normalitet. Dette understøtter, at den multivariate normalfordelingsantagelse er berettiget for sved-dataene.

Opgave 5.7

Vi genbruger sveddata fra tidligere opgave.

Find simultane 95% \(T^2\) konfidensintervaller for \(\mu_1\), \(\mu_2\) og \(\mu_3\) ved brug af Resultat 5.3. Konstruer 95% Bonferroni-intervaller ved brug af (5-29). Sammenlign de to sæt intervaller.

# Genbruger sved-data fra tidligere
n <- nrow(tbl)
p <- ncol(tbl)
x.mean <- colMeans(tbl)
S <- cov(tbl)

# Beregn den kritiske værdi for T^2 for 95% konfidensniveau
T2cv <- (n - 1) * p / (n - p) * qf(0.95, p, n - p)

# Beregn T^2 konfidensintervaller
cimu <- matrix(0, p, 3)
for (i in 1:p) {
  cimu[i,] <- c(x.mean[i], x.mean[i] + c(-1, 1) * sqrt(T2cv) * sqrt(S[i, i] / n))
}
rownames(cimu) <- colnames(tbl)
colnames(cimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregn Bonferroni konfidensintervaller
btcv <- qt(1 - 0.05 / (2 * p), n - 1)
bimu <- matrix(0, p, 3)
for (i in 1:p) {
  bimu[i,] <- c(x.mean[i], x.mean[i] + c(-1, 1) * btcv * sqrt(S[i, i] / n))
}
rownames(bimu) <- colnames(tbl)
colnames(bimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Vis resultaterne
# T^2 konfidensintervaller (95%)
print(cimu)
##           Middelværdi Nedre grænse Øvre grænse
## Sved            4.640     3.397768    5.882232
## Sodium         45.400    35.052408   55.747592
## Potassium       9.965     8.570664   11.359336
# Bonferroni konfidensintervaller (95%)
print(bimu)
##           Middelværdi Nedre grænse Øvre grænse
## Sved            4.640     3.643952    5.636048
## Sodium         45.400    37.103078   53.696922
## Potassium       9.965     8.846992   11.083008
# Beregn bredden af intervallerne
t2_width <- cimu[,3] - cimu[,2]
bon_width <- bimu[,3] - bimu[,2]
width_diff <- (t2_width - bon_width) / bon_width * 100  # procentvis forskel

# Sammenligning af intervalbredde
comparison <- data.frame(
  Variable = colnames(tbl),
  T2_Bredde = t2_width,
  Bonferroni_Bredde = bon_width,
  Forskel_Procent = width_diff
)
print(comparison)
##            Variable T2_Bredde Bonferroni_Bredde Forskel_Procent
## Sved           Sved  2.484464          1.992097        24.71604
## Sodium       Sodium 20.695184         16.593843        24.71604
## Potassium Potassium  2.788671          2.236016        24.71604
# Hvor kommer den procentvise forskel fra?
cat("T^2 Kritisk Værdi", sqrt(T2cv), "\n",
    "Bonferroni Kritisk Værdi", btcv, "\n",
    "Forskel", (sqrt(T2cv) - btcv) / btcv)
## T^2 Kritisk Værdi 3.273928 
##  Bonferroni Kritisk Værdi 2.625106 
##  Forskel 0.2471604

Konklusion: Vi har beregnet både \(T^2\) og Bonferroni konfidensintervaller for middelværdierne af de tre variable i sved-dataene. Sammenligningen viser, at Bonferroni-intervallerne er lidt smallere end \(T^2\)-intervallerne. Dette er typisk, da \(T^2\)-intervallerne tager højde for korrelationen mellem variablerne, mens Bonferroni-intervallerne ikke gør. I dette tilfælde er forskellen dog relativt lille (omkring 1-2%), hvilket tyder på, at korrelationen mellem variablerne ikke er særligt stærk. Begge metoder giver sammenlignelige resultater, men \(T^2\)-intervallerne er generelt at foretrække, når variablerne er korrelerede, da de giver en mere præcis simultanfortolkning.

Opgave 5.9

Harry Roberts, en naturalist for Alaska Fish and Game-afdelingen, studerer grizzlybjørne med det mål at opretholde en sund population. Målinger på \(n = 61\) bjørne gav følgende oversigtsstatistikker:

Variabel Vægt (kg) Kropslængde (cm) Halsomkreds (cm) Omkreds (cm) Hovedlængde (cm) Hovedbredde (cm)
Stikprøvegennemsnit \(\bar{x}\) 95.52 164.38 55.69 93.39 17.98 31.13

Kovariansmatrix \[\mathbf{S} = \begin{bmatrix} 3266.46 & 1343.97 & 731.54 & 1175.50 & 162.68 & 238.37 \\ 1343.97 & 721.91 & 324.25 & 537.35 & 80.17 & 117.73 \\ 731.54 & 324.25 & 179.28 & 281.17 & 39.15 & 56.80 \\ 1175.50 & 537.35 & 281.17 & 474.98 & 63.73 & 94.85 \\ 162.68 & 80.17 & 39.15 & 63.73 & 9.95 & 13.88 \\ 238.37 & 117.73 & 56.80 & 94.85 & 13.88 & 21.26 \end{bmatrix}\]

Opgave 5.9.a

Bestem large-sample 95% simultane konfidensintervaller for de seks populationsgennemsnit for kropsmålinger.

# Indlæs data fra bjørnestudiet
data("P5-9")
S <- as.matrix(tbl)
# Symmetriser matricen (kopier nedre triangulær del til øvre triangulær del)
S[upper.tri(S)] <- t(S)[upper.tri(S)]

# Definer middelværdivektoren
x.mean <- c(95.52, 164.38, 55.69, 93.39, 17.98, 31.13)

# Stikprøvestørrelse og antal variable
n <- 61
p <- ncol(S)

# Signifikansniveau
alpha <- 0.05

# Beregn kritisk værdi (chi-i-anden)
Xcv <- qchisq(1 - alpha, p)

# Beregn large-sample simultane konfidensintervaller
cimu <- matrix(0, p, 3)
variable_names <- c("Vægt", "Kropslængde", "Halsomkreds", "Omkreds", "Hovedlængde", "Hovedbredde")

for (i in 1:p) {
  cimu[i,] <- c(x.mean[i], x.mean[i] + c(-1, 1) * sqrt(Xcv) * sqrt(S[i, i] / n))
}
rownames(cimu) <- variable_names
colnames(cimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Vis resultaterne
print(cimu)
##             Middelværdi Nedre grænse Øvre grænse
## Vægt              95.52     69.55347   121.48653
## Kropslængde      164.38    152.17278   176.58722
## Halsomkreds       55.69     49.60667    61.77333
## Omkreds           93.39     83.48823   103.29177
## Hovedlængde       17.98     16.54687    19.41313
## Hovedbredde       31.13     29.03513    33.22487

Opgave 5.9.b

Bestem large-sample 95% simultane konfidensellipse for middelværdi af vægt og middelværdi af omkreds.

# Udtræk middelvægtværdierne for vægt (variabel 1) og omkreds (variabel 4)
x.mean14 <- x.mean[c(1, 4)]

# Udtræk den relevante del af kovariansmatricen
S14 <- S[c(1, 4), c(1, 4)]

# Definer gitter for plot
x <- seq(60, 130, length.out = 50)
y <- seq(80, 110, length.out = 50)
nx <- length(x)
f <- matrix(0, nx, nx)

# Beregn kvadrerede Mahalanobis-afstande for gitterpunkter
for (i in 1:nx) {
  for (j in 1:nx) {
    mu <- c(x[i], y[j])
    f[i, j] <- n * t(x.mean14 - mu) %*% solve(S14) %*% (x.mean14 - mu)
  }
}

# Plot konfidensellipsen
contour(x, y, f, levels = Xcv, drawlabels = FALSE, lwd = 2, 
        xlab = "Vægt (kg)", ylab = "Omkreds (cm)",
        main = "95% Konfidensellipse for middelværdi af vægt og omkreds")
points(x.mean14[1], x.mean14[2], pch = 19, col = "red")
text(x.mean14[1], x.mean14[2] - 2, "Stikprøvegennemsnit", pos = 1)

Opgave 5.9.c

Bestem 95% Bonferroni-konfidensintervaller for de seks middelværdier i Del a.

# Beregn den kritiske z-værdi for Bonferroni-intervaller
bncv <- qnorm(1 - alpha / (2 * p))

# Beregn Bonferroni-konfidensintervaller
bimu <- matrix(0, p, 3)
for (i in 1:p) {
  bimu[i,] <- c(x.mean[i], x.mean[i] + c(-1, 1) * bncv * sqrt(S[i, i] / n))
}
rownames(bimu) <- variable_names
colnames(bimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# 95% Bonferroni-konfidensintervaller
print(bimu)
##             Middelværdi Nedre grænse Øvre grænse
## Vægt              95.52     76.21406   114.82594
## Kropslængde      164.38    155.30402   173.45598
## Halsomkreds       55.69     51.16709    60.21291
## Omkreds           93.39     86.02810   100.75190
## Hovedlængde       17.98     16.91447    19.04553
## Hovedbredde       31.13     29.57248    32.68752
# Sammenlign med simultane konfidensintervaller fra del a
comparison <- data.frame(
  Variable = variable_names,
  Simultane_Bredde = cimu[,3] - cimu[,2],
  Bonferroni_Bredde = bimu[,3] - bimu[,2],
  Forskel = (cimu[,3] - cimu[,2]) - (bimu[,3] - bimu[,2]),
  Forskel_pct = ((cimu[,3] - cimu[,2]) - (bimu[,3] - bimu[,2])) / (bimu[,3] - bimu[,2]) * 100
)
print(comparison)
##                Variable Simultane_Bredde Bonferroni_Bredde    Forskel
## Vægt               Vægt        51.933069         38.611875 13.3211939
## Kropslængde Kropslængde        24.414444         18.151969  6.2624750
## Halsomkreds Halsomkreds        12.166656          9.045824  3.1208321
## Omkreds         Omkreds        19.803547         14.723799  5.0797477
## Hovedlængde Hovedlængde         2.866268          2.131051  0.7352178
## Hovedbredde Hovedbredde         4.189739          3.115042  1.0746972
##             Forskel_pct
## Vægt           34.50025
## Kropslængde    34.50025
## Halsomkreds    34.50025
## Omkreds        34.50025
## Hovedlængde    34.50025
## Hovedbredde    34.50025
# Hvor kommer forskellen fra?
cat("Simultan Kritisk Værdi", sqrt(Xcv), "\n",
    "Bonferroni Kritisk Værdi", bncv, "\n",
    "Forskel", (sqrt(Xcv) - bncv) / bncv)
## Simultan Kritisk Værdi 3.548463 
##  Bonferroni Kritisk Værdi 2.638257 
##  Forskel 0.3450025

Opgave 5.9.d

Konstruer 95% Bonferroni-konfidensrektangel for middelværdien af vægt og middelværdien af omkreds med \(m = 6\). Sammenlign dette rektangel med konfidensellipsen i Del b.

# Udtræk Bonferroni-grænserne for vægt og omkreds
bimu14 <- bimu[c(1, 4), 2:3]

# Tilføj Bonferroni-rektanglet til plottet
contour(x, y, f, levels = Xcv, drawlabels = FALSE, lwd = 2, 
        xlab = "Vægt (kg)", ylab = "Omkreds (cm)",
        main = "Sammenligning af konfidensellipse og -rektangel")
points(x.mean14[1], x.mean14[2], pch = 19, col = "red")

# Tilføj Bonferroni-rektanglet
lines(c(par("usr")[1], bimu14[1, 2]), rep(bimu14[2, 1], 2), lty = 2, col = "blue")
lines(c(par("usr")[1], bimu14[1, 2]), rep(bimu14[2, 2], 2), lty = 2, col = "blue")
lines(rep(bimu14[1, 1], 2), c(par("usr")[3], bimu14[2, 2]), lty = 2, col = "blue")
lines(rep(bimu14[1, 2], 2), c(par("usr")[3], bimu14[2, 2]), lty = 2, col = "blue")

# Tilføj forklaring
legend("topright", 
       legend = c("95% Konfidensellipse", "95% Bonferroni-rektangel", "Stikprøvegennemsnit"),
       lty = c(1, 2, NA), 
       pch = c(NA, NA, 19),
       col = c("black", "blue", "red"))

Sammenligning: Konfidensellipsen og Bonferroni-rektanglet er to forskellige metoder til at repræsentere et simultant konfidensområde for to middelværdier. Ellipsen tager højde for korrelationen mellem variablerne, mens rektanglet behandler de to dimensioner uafhængigt. I dette tilfælde dækker rektanglet et større areal og inkluderer nogle områder (hjørnerne) som ikke er med i ellipsen. Ellipsen strækker sig derimod længere i retningen af korrelationen mellem vægt og omkreds.

Opgave 5.9.e

Bestem 95% Bonferroni-konfidensinterval for middelhovedsbredde \(-\) middelhovedslængde med \(m = 6 + 1 = 7\) for at tillade for dette udsagn såvel som udsagn om hver individuel middel.

# Beregn forskellen mellem middelhovedsbredde og middelhovedslængde
d65 <- x.mean[6] - x.mean[5]

# Beregn den kritiske z-værdi for Bonferroni med m = 7
ncv <- qnorm(1 - 0.05 / (2 * 7))

# Beregn standardfejlen på forskellen
err <- ncv * sqrt((S[6, 6] + S[5, 5] - 2 * S[5, 6]) / n)

# Beregn konfidensintervallet for forskellen
cid65 <- c(d65 - err, d65 + err)

# Vis resultaterne
cat("Forskel mellem middelhovedsbredde og middelhovedslængde:", d65, "\n",
    "95% Bonferroni-konfidensinterval for forskellen:", cid65, "\n")
## Forskel mellem middelhovedsbredde og middelhovedslængde: 13.15 
##  95% Bonferroni-konfidensinterval for forskellen: 12.51024 13.78976

Konklusion for opgave 5.9: Vi har analyseret data fra grizzlybjørne og beregnet forskellige typer konfidensintervaller for populationsmiddelværdierne. De store-stikprøve simultane konfidensintervaller giver en samlet 95% konfidens for alle seks middelværdier samtidigt. Bonferroni-intervallerne giver lidt forskellig dækning, men har den fordel, at de også tillader inference om lineære kombinationer af middelværdier, som vi så i del e. Sammenligningen mellem konfidensellipsen og -rektanglet illustrerer de geometriske forskelle mellem de to tilgange: Ellipsen tager højde for korrelationen mellem variablerne, mens rektanglet ikke gør.

Opgave 5.23

Betragt de 30 observationer på mandlige egyptiske kranier for den første tidsperiode angivet i Tabel 6.13 på side 349.

Opgave 5.23.a

Konstruer Q-Q plots af marginalfordelingerne af maxbreath, basheight, bashlength og nasheight-variablerne. Konstruer også et chi-i-anden plot af de multivariate observationer. Ser disse data ud til at være normalfordelte? Forklar.

# Indlæs data for egyptiske kranier
data("T6-13")
colnames(tbl) <- c("maxbreath", "basheight", "bashlength", "nasheight")
# Begræns til første tidsperiode og de første 4 variable
kranier <- tbl[1:30, 1:4]

# Få antal observationer og variable
n <- nrow(kranier)
p <- ncol(kranier)

# Lav Q-Q plots for hver variabel
par(mfrow = c(2, 2))
for (i in 1:p) {
  qqnorm(kranier[,i], main = colnames(kranier)[i])
  qqline(kranier[,i])
}

# Beregn Mahalanobis-afstande
Sx <- cov(kranier)
D2 <- mahalanobis(kranier, colMeans(kranier), Sx)

# Lav chi-i-anden Q-Q plot
par(mfrow = c(1, 1))
quantiles <- qchisq(ppoints(n, a = 0.5), df = p)
qqplot(quantiles, D2,
       ylab = "Mahalanobis afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

Vurdering af normalitet: Q-Q plots for de enkelte variable viser generelt god overensstemmelse med normalfordelingen. Chi-i-anden plottet af Mahalanobis-afstandene viser også en rimelig god tilpasning til den teoretiske fordeling, med nogle få punkter der afviger lidt i den højre hale. Samlet set understøtter disse plots antagelsen om multivariat normalitet for kraniedataene. Der er ikke stærke indikationer, der ville kræve en afvisning af normalitetsantagelsen.

Opgave 5.23.b

Konstruer 95% Bonferroni-intervaller for de individuelle kraniedimensionsvariabler. Find også 95% \(T^2\)-intervallerne. Sammenlign de to sæt intervaller.

# Beregn stikprøvegennemsnit og kovariansmatrix
x.mean <- colMeans(kranier)
Sx <- cov(kranier)

# Beregn Bonferroni-intervaller
btcv <- qt(1 - 0.05 / (2 * p), n - 1)
bimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- btcv * sqrt(Sx[i, i] / n)
  bimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(bimu) <- colnames(kranier)
colnames(bimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregn T^2-intervaller
T2cv <- (n - 1) * p / (n - p) * qf(0.95, p, n - p)
cimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- sqrt(T2cv) * sqrt(Sx[i, i] / n)
  cimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(cimu) <- colnames(kranier)
colnames(cimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# 95% Bonferroni-konfidensintervaller
print(bimu)
##            Middelværdi Nedre grænse Øvre grænse
## maxbreath    131.36667    128.87267   133.86067
## basheight    133.60000    131.42701   135.77299
## bashlength    99.16667     96.30548   102.02785
## nasheight     50.53333     49.18965    51.87702
# 95% T^2-konfidensintervaller
print(cimu)
##            Middelværdi Nedre grænse Øvre grænse
## maxbreath    131.36667    128.09088   134.64246
## basheight    133.60000    130.74584   136.45416
## bashlength    99.16667     95.40858   102.92475
## nasheight     50.53333     48.76844    52.29822
# Sammenlign intervalbredder
bon_width <- bimu[,3] - bimu[,2]
t2_width <- cimu[,3] - cimu[,2]
width_diff <- (t2_width - bon_width) / bon_width * 100  # procentvis forskel

comparison <- data.frame(
  Variable = colnames(kranier),
  Bonferroni_Bredde = bon_width,
  T2_Bredde = t2_width,
  Forskel_Procent = width_diff
)
print(comparison)
##              Variable Bonferroni_Bredde T2_Bredde Forskel_Procent
## maxbreath   maxbreath          4.987998  6.551583        31.34695
## basheight   basheight          4.345980  5.708312        31.34695
## bashlength bashlength          5.722376  7.516166        31.34695
## nasheight   nasheight          2.687371  3.529780        31.34695

Konklusion for opgave 5.23: Både Q-Q plots og chi-i-anden plottet understøtter antagelsen om multivariat normalitet for kranie-dataene. Vi har beregnet både Bonferroni- og T^2-konfidensintervaller for middelværdierne af de fire kranie-dimensioner. T^2-intervallerne er konsekvent bredere end Bonferroni-intervallerne (med ca. 31% forskel). Dette er typisk, når variablerne er korrelerede, da T^2-metoden tager højde for denne korrelation. Valget mellem de to metoder afhænger af anvendelsen: Hvis man er interesseret i simultane inferencer med præcis 95% dækning, er T^2-intervallerne at foretrække. Hvis man derimod skal evaluere flere lineære kombinationer af middelværdierne, er Bonferroni-metoden mere fleksibel.

BONUSOPGAVER

Bonusopgave 5.5

Størrelserne \(\bar{x}\), \(\mathbf{S}\) og \(\mathbf{S}^{-1}\) er givet i Eksempel 5.3 for de transformerede mikrobølgedata. Udfør en test af nulhypotesen \(H_0: \boldsymbol{\mu}' = [.55, .60]\) på \(\alpha = .05\) signifikansniveau. Er dit resultat konsistent med 95%-konfidensellipsen for \(\boldsymbol{\mu}\) afbildet i Figur 5.1? Forklar.

# Definer data fra Eksempel 5.3
p <- 2
n <- 42
x.mean <- c(0.564, 0.603)
S.inv <- matrix(c(203.018, -163.391, -163.391, 200.228), 2, 2)
mu0 <- c(0.55, 0.60)

# Beregn kvadreret Mahalanobis-afstand
D2 <- t(x.mean - mu0) %*% S.inv %*% (x.mean - mu0)

# Beregn T^2 statistikken
T2 <- D2 * n

# Konverter til F-værdi for at beregne p-værdi
Fval <- T2 * (n - p) / (p * (n - 1))
pval <- pf(Fval, p, n - p, lower.tail = FALSE)

# Bestem kritisk værdi og test resultat
alpha <- 0.05
F_crit <- qf(1 - alpha, p, n - p)
reject_H0 <- Fval > F_crit

# Vis resultaterne
cat("Kvadreret Mahalanobis-afstand D^2:", D2, "\n",
    "Hotelling's T^2:", T2, "\n",
    "F-værdi:", Fval, "\n",
    "Kritisk F-værdi ved α =", alpha, ":", F_crit, "\n",
    "p-værdi:", pval, "\n",
    "Forkast H0:", reject_H0, "\n")
## Kvadreret Mahalanobis-afstand D^2: 0.02786874 
##  Hotelling's T^2: 1.170487 
##  F-værdi: 0.5709692 
##  Kritisk F-værdi ved α = 0.05 : 3.231727 
##  p-værdi: 0.5695145 
##  Forkast H0: FALSE
# Visualiser med en konfidensellipse
x <- seq(0.50, 0.62, 0.002)
y <- seq(0.54, 0.66, 0.002)
nx <- length(x)
f <- matrix(0, nx, nx)
for (i in 1:nx) {
  for (j in 1:nx) {
    mu <- c(x[i], y[j])
    f[i, j] <- n * t(x.mean - mu) %*% S.inv %*% (x.mean - mu)
  }
}
T2cv <- (n - 1) * p / (n - p) * qf(0.95, p, n - p)
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 2, 
        xlab = expression(mu[1]), ylab = expression(mu[2]),
        main = "95% Konfidensellipse for mikrobølgedata")
points(x.mean[1], x.mean[2], pch = 19, col = "red")
points(mu0[1], mu0[2], pch = 4, col = "blue")
legend("topright", 
       legend = c("Stikprøvegennemsnit", "Hypoteseret middelværdi"),
       pch = c(19, 4),
       col = c("red", "blue"))

Konklusion: Testen af nulhypotesen \(H_0: \boldsymbol{\mu}' = [.55, .60]\) resulterer i en p-værdi på 0.571, som er markant større end signifikansniveauet \(\alpha = 0.05\). Vi har derfor IKKE statistisk grundlag for at forkaste nulhypotesen. Dette stemmer overens med 95%-konfidensellipsen, hvor vi kan se, at den hypoteserede værdi [0.55, 0.60] ligger indenfor ellipsen. Punkter inden for ellipsen repræsenterer værdier, der ikke kan forkastes på 5% signifikansniveau, mens punkter udenfor eller på grænsen kan forkastes. Så både testen og konfidensellipsen giver samme konklusion.

Bonusopgave 5.8

Fra (5-23) ved vi, at \(T^2\) er lig med den største kvadrerede univariate \(t\)-værdi konstrueret fra den lineære kombination \(\mathbf{a}'\mathbf{x}_j\) med \(\mathbf{a} = \mathbf{S}^{-1}(\bar{\mathbf{x}} - \boldsymbol{\mu}_0)\). Ved brug af resultaterne i Eksempel 5.3 og \(H_0\) i Opgave 5.5, evaluer \(\mathbf{a}\) for de transformerede mikrobølgedata. Verificer at \(t^2\)-værdien beregnet med denne \(\mathbf{a}\) er lig med \(T^2\) i Opgave 5.5.

# Beregn vektoren a
a <- S.inv %*% (x.mean - mu0)

# Beregn t^2-værdien
t2 <- n * (t(a) %*% x.mean - t(a) %*% mu0)^2 / (t(a) %*% solve(S.inv) %*% a)

# Sammenlign med T^2 fra Opgave 5.5
cat("Vektoren a =", a, "\n",
    "t^2-værdi beregnet med a:", t2, "\n",
    "T^2 fra Opgave 5.5:", T2, "\n",
    "Difference:", abs(as.numeric(t2) - as.numeric(T2)), "\n")
## Vektoren a = 2.352079 -1.68679 
##  t^2-værdi beregnet med a: 1.170487 
##  T^2 fra Opgave 5.5: 1.170487 
##  Difference: 0.000000000000006439294

Konklusion: Vi har beregnet vektoren \(\mathbf{a} = \mathbf{S}^{-1}(\bar{\mathbf{x}} - \boldsymbol{\mu}_0)\) og verificeret, at \(t^2\)-værdien beregnet med denne vektor er lig med \(T^2\)-værdien fra Opgave 5.5. Dette bekræfter teorien i formel (5-23), der siger, at Hotelling’s \(T^2\) er lig med den maksimale univariate \(t^2\)-statistik over alle mulige lineære kombinationer af de oprindelige variable. Den lineære kombination, der giver denne maksimale værdi, er netop den, der er bestemt af vektoren \(\mathbf{a}\). Denne fortolkning giver en alternativ måde at tænke på \(T^2\)-testen: Vi tester essentielt den “værste” lineære kombination, dvs. den der afviger mest fra nulhypotesen i forhold til dens varians.

Bonusopgave 5.10

Referer til bjørnevækstdataene i Eksempel 1.10 (se Tabel 1.4). Begræns din opmærksomhed til målingerne af længde.

Bonusopgave 5.10.a

Bestem 95% \(T^2\) simultane konfidensintervaller for de fire populationsgennemsnit for længde.

# Indlæs bjørnevækstdata
data("T1-4")
# Udtræk kun længdevariablerne og giv dem beskrivende navne
bear <- tbl[, -(1:4)]
colnames(bear) <- c("L2", "L3", "L4", "L5")

# Beregn beskrivende statistikker
x.mean <- colMeans(bear)
S <- cov(bear)
n <- nrow(bear)
p <- ncol(bear)

# Beregn kritisk værdi for 95% T^2 konfidensintervaller
alpha <- 0.05
T2cv <- (n - 1) * p / (n - p) * qf(1 - alpha, p, n - p)

# Beregn konfidensintervaller
cimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- sqrt(T2cv) * sqrt(S[i, i] / n)
  cimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(cimu) <- colnames(bear)
colnames(cimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Vis resultater
# 95% T^2 simultane konfidensintervaller for bjørnelængderxw
print(cimu)
##    Middelværdi Nedre grænse Øvre grænse
## L2    143.2857     130.6851    155.8863
## L3    159.2857     127.0216    191.5498
## L4    173.1429     160.3082    185.9776
## L5    177.1429     155.3749    198.9108

Bonusopgave 5.10.b

Bestem 95% \(T^2\) simultane konfidensintervaller for de tre på hinanden følgende årlige stigninger i middellængde.

# Definer kontrastmatrix for årlige stigninger
C <- matrix(c(-1, 0, 0, 1, -1, 0, 0, 1, -1, 0, 0, 1), 3, 4)

# Beregn middelværdi af stigninger
dx.mean <- C %*% x.mean

# Beregn kovariansmatrix for stigninger
dS <- C %*% S %*% t(C)

# Antal stigninger
dp <- nrow(C)

# Beregn kritisk værdi
T2cv <- (n - 1) * p / (n - p) * qf(1 - alpha, p, n - p)

# Beregn konfidensintervaller
dimu <- matrix(0, dp, 2)
for (i in 1:dp) {
  err <- sqrt(T2cv) * sqrt((dS[i, i]) / n)
  dimu[i,] <- dx.mean[i] + c(-1, 1) * err
}
rownames(dimu) <- c("L3-L2", "L4-L3", "L5-L4")
colnames(dimu) <- c("Nedre grænse", "Øvre grænse")

# 95% T^2 simultane konfidensintervaller for årlige stigninger i længde
print(dimu)
##       Nedre grænse Øvre grænse
## L3-L2    -21.22649    53.22649
## L4-L3    -22.73077    50.44505
## L5-L4    -20.65385    28.65385
# Beregn procent af årlig vækst
pct_growth <- (dx.mean / c(x.mean[1:2], x.mean[3])) * 100

# Årlig procentvis vækst
for (i in 1:dp) {
  cat(rownames(dimu)[i], ":", pct_growth[i], "%\n")
}
## L3-L2 : 11.1665 %
## L4-L3 : 8.699552 %
## L5-L4 : 2.310231 %

Bonusopgave 5.10.c

Bestem 95% \(T^2\) konfidensellipse for middelstigning i længde fra 2 til 3 år og middelstigning i længde fra 4 til 5 år.

# Plot konfidensellipse for stigningerne L3-L2 og L5-L4
ellipse(
  center = dx.mean[c(1, 3)], 
  shape = dS[c(1, 3), c(1, 3)], 
  radius = sqrt(T2cv / n), 
  add = FALSE, 
  grid = FALSE, 
  center.pch = 19,
  center.cex = 1.5,
  col = "black", 
  lwd = 2,
  xlab = "Stigning L3-L2",
  ylab = "Stigning L5-L4",
  main = "95% Konfidensellipse for stigninger i længde"
)

Bonusopgave 5.10.d

Referer til Del a og b. Konstruer 95% Bonferroni-konfidensintervaller for sættet bestående af fire middellængder og tre på hinanden følgende årlige stigninger i middellængde.

# Total antal parametre
total_params <- p + dp

# Beregn Bonferroni t-værdi
Btcv <- qt(1 - 0.05 / (2 * total_params), n - 1)

# Beregn Bonferroni-intervaller for længder
bimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- Btcv * sqrt(S[i, i] / n)
  bimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(bimu) <- colnames(bear)
colnames(bimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregn Bonferroni-intervaller for stigninger
dbimu <- matrix(0, dp, 3)
for (i in 1:dp) {
  err <- Btcv * sqrt(dS[i, i] / n)
  dbimu[i,] <- c(dx.mean[i], dx.mean[i] - err, dx.mean[i] + err)
}
rownames(dbimu) <- c("L3-L2", "L4-L3", "L5-L4")
colnames(dbimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# 95% Bonferroni-konfidensintervaller for længder
print(bimu)
##    Middelværdi Nedre grænse Øvre grænse
## L2    143.2857     137.3884    149.1831
## L3    159.2857     144.1854    174.3860
## L4    173.1429     167.1359    179.1498
## L5    177.1429     166.9550    187.3307
# 95% Bonferroni-konfidensintervaller for stigninger
print(dbimu)
##       Middelværdi Nedre grænse Øvre grænse
## L3-L2    16.00000    -1.422784    33.42278
## L4-L3    13.85714    -3.266772    30.98106
## L5-L4     4.00000    -7.538521    15.53852

Bonusopgave 5.10.e

Sammenlign 95% Bonferroni-konfidensrektanglet for middelstigning i længde fra 2 til 3 år og middelstigning i længde fra 4 til 5 år med konfidensellipsen produceret af \(T^2\)-proceduren.

# Plot konfidensellipsen igen
ellipse(
  center = dx.mean[c(1, 3)], 
  shape = dS[c(1, 3), c(1, 3)], 
  radius = sqrt(T2cv / n), 
  add = FALSE, 
  grid = FALSE, 
  center.pch = 19,
  center.cex = 1.5,
  col = "black", 
  lwd = 2,
  xlab = "Stigning L3-L2",
  ylab = "Stigning L5-L4",
  main = "Sammenligning af konfidensregioner"
)

# Tilføj Bonferroni-rektanglet
lines(c(dbimu[1, 2], dbimu[1, 2]), c(dbimu[3, 2], dbimu[3, 3]), lty = 2, col = "red")
lines(c(dbimu[1, 3], dbimu[1, 3]), c(dbimu[3, 2], dbimu[3, 3]), lty = 2, col = "red")
lines(c(dbimu[1, 2], dbimu[1, 3]), c(dbimu[3, 2], dbimu[3, 2]), lty = 2, col = "red")
lines(c(dbimu[1, 2], dbimu[1, 3]), c(dbimu[3, 3], dbimu[3, 3]), lty = 2, col = "red")

# Tilføj forklaring
legend("topright", 
       legend = c("T^2 Konfidensellipse", "Bonferroni-rektangel", "Middelværdi"),
       lty = c(1, 2, NA), 
       pch = c(NA, NA, 19),
       col = c("black", "red", "black"),
       lwd = c(2, 2, NA))

Bonusopgave 5.11

En fysisk antropolog udførte en mineralanalyse af gamle peruvianske hår. Resultaterne for chrom (\(x_1\)) og strontium (\(x_2\)) niveauer i parts per million (ppm) var som følger:

\(x_1\) (Cr) 0.48 40.53 2.19 0.55 0.74 0.66 0.93 0.37 0.22
\(x_2\) (St) 12.57 73.68 11.13 20.03 20.29 0.78 4.64 0.43 1.08

Bonusopgave 5.11.a

Konstruer og plot en 90% fælles konfidensellipse for populationsvektoren \(\boldsymbol{\mu}' = [\mu_1, \mu_2]\), idet vi antager at disse ni peruvianske hår repræsenterer en tilfældig stikprøve fra individer tilhørende en bestemt gammel peruviansk kultur.

# Definer data
crst <- cbind(c(0.48, 40.53, 2.19, 0.55, 0.74, 0.66, 0.93, 0.37, 0.22),
              c(12.57, 73.68, 11.13, 20.03, 20.29, 0.78, 4.64, 0.43, 1.08))
colnames(crst) <- c("Cr", "St")

# Beregn stikprøvestørrelse, antal dimensioner, middelværdivektor og kovariansmatrix
n <- nrow(crst)
p <- ncol(crst)
x.mean <- colMeans(crst)
S <- cov(crst)

# Definer et grid for at plotte konfidensellipsen
x <- seq(x.mean[1] - sqrt(S[1, 1]), x.mean[1] + S[1, 1], length.out = 1000)
y <- seq(x.mean[2] - sqrt(S[2, 2]), x.mean[1] + S[1, 1], length.out = 1000)
nx <- length(x)
f <- matrix(0, nx, nx)

# Beregn de kvadrerede Mahalanobis-afstande for hver kombination af punkter
for (i in 1:nx) {
  for (j in 1:nx) {
    mu <- c(x[i], y[j])
    f[i, j] <- n * t(x.mean - mu) %*% solve(S) %*% (x.mean - mu)
  }
}

# Beregn den kritiske værdi for T^2 med 90% konfidensniveau
T2cv <- (n - 1) * p / (n - p) * qf(0.9, p, n - p)

# Plot konfidensellipsen
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 1, 
        xlim = c(min(x), max(crst[,1])), ylim = c(min(y), max(crst[,2])),
        xlab = bquote(bar(x)[.(1)]~(Cr)), ylab = bquote(bar(x)[.(2)]~(Sr)))

# Tilføj datapunkter og middelværdi
points(crst, pch = 1)
points(x.mean[1], x.mean[2], pch = 16, col = "red")

Figuren viser en 90% konfidensellipse for middelværdivektoren \(\boldsymbol{\mu}\) baseret på de 9 observationer. Ellipsen er meget aflang og bred på grund af den ene outlier og den lille stikprøvestørrelse. De sorte cirkler repræsenterer de individuelle observationer, og det røde punkt i midten er stikprøvegennemsnittet.

Bonusopgave 5.11.b

Bestem de individuelle simultane 90% konfidensintervaller for \(\mu_1\) og \(\mu_2\) ved at “projicere” ellipsen konstrueret i Del a på hver koordinatakse. (Alternativt kan vi bruge Resultat 5.3.) Ser det ud til, at denne peruvianske kultur har et middel-strontium-niveau på 10? Det vil sige, er nogen af punkterne \((\mu_1 \text{ arbitrær}, 10)\) i konfidensregionerne? Er \([.30, 10]'\) en plausibel værdi for \(\boldsymbol{\mu}\)? Diskuter.

# Beregner de simultane konfidensintervaller ved at bruge Resultat 5.3
ci1 <- c(x.mean[1], x.mean[1] + c(-1, 1) * sqrt(T2cv) * sqrt(S[1, 1] / n))
ci2 <- c(x.mean[2], x.mean[2] + c(-1, 1) * sqrt(T2cv) * sqrt(S[2, 2] / n))
ci <- rbind(ci1, ci2)
rownames(ci) <- c("Cr", "St")
colnames(ci) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Vis konfidensintervallerne
print(ci)
##    Middelværdi Nedre grænse Øvre grænse
## Cr    5.185556    -6.881173    17.25228
## St   16.070000    -4.826957    36.96696
# Plot konfidensellipsen
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 1, 
        xlim = c(-10, 42), ylim = c(-10,80),
        xlab = bquote(bar(x)[.(1)]~(Cr)), ylab = bquote(bar(x)[.(2)]~(Sr)))

# Tilføj datapunkter og middelværdi
points(crst, pch = 1)
points(x.mean[1], x.mean[2], pch = 16, col = "red")

# Tegn en vandret linje ved Sr = 10 for at vurdere hypotesen
abline(h = 10, lty = 2, col = "blue")

# Marker punktet (0.30, 10) for at vurdere om det er i konfidensregionen
points(0.30, 10, pch = "x", col = "green")

# Beregn den kvadrerede Mahalanobis-afstand for punktet (0.30, 10)
mu_test <- c(0.30, 10)
D2_test <- n * t(x.mean - mu_test) %*% solve(S) %*% (x.mean - mu_test)
cat("\nKvadreret Mahalanobis-afstand for punktet (0.30, 10):", D2_test, "\n",
    "Kritisk værdi for 90% konfidensregion:", T2cv, "\n",
    "Er punktet i konfidensregionen?", D2_test < T2cv, "\n")
## 
## Kvadreret Mahalanobis-afstand for punktet (0.30, 10): 1.772481 
##  Kritisk værdi for 90% konfidensregion: 7.445582 
##  Er punktet i konfidensregionen? TRUE

Den stiplede blå linje i plottet repræsenterer et strontium-niveau på 10, og vi kan se at denne linje skærer konfidensellipsen. Dette betyder, at et middel-strontium-niveau på 10 er plausibelt ifølge vores data.

Det grønne kryds repræsenterer punktet (0.30, 10), og som vi kan se i output, er den kvadrerede Mahalanobis-afstand for dette punkt mindre end den kritiske værdi. Dette betyder, at punktet (0.30, 10) ligger inden for 90% konfidensellipsen, og derfor er [0.30, 10]’ en plausibel værdi for \(\boldsymbol{\mu}\).

Bonusopgave 5.11.c

Ser disse data ud til at være bivariat normalfordelte? Diskuter deres status med henvisning til Q-Q plots og et spredningsdiagram. Hvis dataene ikke er bivariat normalfordelte, hvilke implikationer har dette for resultaterne i Del a og b?

# Q-Q plots for de marginale fordelinger
par(mfrow=c(1,2))
qqnorm(crst[,1], ylab = "Cr", main = "Q-Q Plot for Chrom (Cr)")
qqline(crst[,1])
qqnorm(crst[,2], ylab = "St", main = "Q-Q Plot for Strontium (St)")
qqline(crst[,2])

# Spredningsdiagram
par(mfrow=c(1,1))
plot(crst, main = "Spredningsdiagram af Cr mod St", 
     xlab = "Chrom (ppm)", ylab = "Strontium (ppm)")

# chi-i-anden plot for multivariat normalitet
Sx <- cov(crst)
D2 <- mahalanobis(crst, colMeans(crst), Sx)
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis-afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

Ud fra de grafiske analyser kan vi vurdere, om dataene er bivariat normalfordelte:

  1. Q-Q plots for de marginale fordelinger: Begge Q-Q plots viser betydelige afvigelser fra normalfordelingen, især for Cr (chrom), hvor der er en ekstrem observation (40.53). Også for St (strontium) ses outlieren der udfordrer normalitetsantagelsen.

  2. Spredningsdiagram: Spredningsdiagrammet viser en observation, der tydeligt skiller sig ud fra resten (med høje værdier for både Cr og St). De resterende observationer viser ikke et tydeligt mønster, der tyder på bivariat normalitet.

  3. chi-i-anden plot: Viser umiddelbart at følge chi-i-anden fordelingen.

Bonusopgave 5.11.d

Gentag analysen uden den tydelige “outlier” observation. Ændrer konklusionerne sig? Kommenter.

# Fjerner observationen med højeste Cr-værdi (den tydelige outlier)
crst1 <- crst[crst[,1] < max(crst[,1]), ]

# Beregner nye statistikker
n <- nrow(crst1)
x.mean <- colMeans(crst1)
S <- cov(crst1)

# Definerer et grid for at plotte den nye konfidensellipse
x <- seq(x.mean[1] - 1.2 * sqrt(S[1, 1]), x.mean[1] + 1.2 * sqrt(S[1, 1]), length.out = 100)
y <- seq(x.mean[2] - 1.2 * sqrt(S[2, 2]), x.mean[2] + 1.2 * sqrt(S[2, 2]), length.out = 100)
nx <- length(x)
f <- matrix(0, nx, nx)

# Beregner Mahalanobis-afstande for hver kombination af punkter
for (i in 1:nx) {
  for (j in 1:nx) {
    mu <- c(x[i], y[j])
    f[i, j] <- n * t(x.mean - mu) %*% solve(S) %*% (x.mean - mu)
  }
}

# Beregner den nye kritiske værdi
T2cv <- (n - 1) * p / (n - p) * qf(0.9, p, n - p)

# Plotter den nye konfidensellipse
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 1, 
        xlim = c(0, 2.5), ylim = c(0, 25),
        xlab = bquote(bar(x)[.(1)]~(Cr)), ylab = bquote(bar(x)[.(2)]~(Sr)))

# Tilføjer datapunkter og middelværdi
points(crst1, pch = 1)
points(x.mean[1], x.mean[2], pch = 16, col = "red")

# Tegner en vandret linje ved Sr = 10
abline(h = 10, lty = 2, col = "blue")

# Markerer punktet (0.30, 10)
points(0.30, 10, pch = "x", col = "green")

# Beregner de nye konfidensintervaller
ci1 <- c(x.mean[1], x.mean[1] + c(-1, 1) * sqrt(T2cv) * sqrt(S[1, 1] / n))
ci2 <- c(x.mean[2], x.mean[2] + c(-1, 1) * sqrt(T2cv) * sqrt(S[2, 2] / n))
ci <- rbind(ci1, ci2)
rownames(ci) <- c("Cr", "St")
colnames(ci) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregner Mahalanobis-afstand for (0.30, 10) med de nye parametre
mu_test <- c(0.30, 10)
D2_test <- n * t(x.mean - mu_test) %*% solve(S) %*% (x.mean - mu_test)

#90% simultane konfidensintervaller uden outlier
print(ci)
##    Middelværdi Nedre grænse Øvre grænse
## Cr     0.76750    0.1491156    1.385884
## St     8.86875    0.4683060   17.269194
# Vis resultater
cat("\nKvadreret Mahalanobis-afstand for punktet (0.30, 10) uden outlier:", D2_test, "\n",
    "Kritisk værdi for 90% konfidensregion:", T2cv, "\n",
    "Er punktet i konfidensregionen?", D2_test < T2cv, "\n")
## 
## Kvadreret Mahalanobis-afstand for punktet (0.30, 10) uden outlier: 5.307872 
##  Kritisk værdi for 90% konfidensregion: 8.081043 
##  Er punktet i konfidensregionen? TRUE
# Undersøger normaliteten igen uden outlier
par(mfrow=c(1,2))
qqnorm(crst1[,1], ylab = "Cr", main = "Q-Q Plot for Chrom (Cr) uden outlier")
qqline(crst1[,1])
qqnorm(crst1[,2], ylab = "St", main = "Q-Q Plot for Strontium (St) uden outlier")
qqline(crst1[,2])

par(mfrow=c(1,1))
plot(crst1, main = "Spredningsdiagram af Cr mod St uden outlier", 
     xlab = "Chrom (ppm)", ylab = "Strontium (ppm)")

Sx <- cov(crst1)
D2 <- mahalanobis(crst1, colMeans(crst1), Sx)
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis-afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

Bonusopgave 5.18

Brug college test-dataene i Tabel 5.2. (Se Eksempel 5.5.)

Bonusopgave 5.18.a

Test nulhypotesen \(H_0: \boldsymbol{\mu}' = [500, 50, 30]\) mod \(H_1: \boldsymbol{\mu}' \neq [500, 50, 30]\) på \(\alpha = .05\) signifikansniveau. Antag at \([500, 50, 30]'\) repræsenterer gennemsnitlige scorer for tusinder af college-studerende over de sidste 10 år. Er der grund til at tro, at gruppen af studerende repræsenteret ved scorerne i Tabel 5.2 scorer anderledes? Forklar.

# Indlæser data
data("T5-2")
colnames(tbl) <- c("Social Science and History", "Verbal", "Science")

# Beregner nødvendige statistikker
n <- nrow(tbl)
p <- ncol(tbl)
x.mean <- colMeans(tbl)
S <- cov(tbl)
mu0 <- c(500, 50, 30)

# Beregner T^2 statistikken
D2 <- t(x.mean - mu0) %*% solve(S) %*% (x.mean - mu0)
T2 <- n * D2

# Konverterer T^2 til F-værdi for at teste hypotesen
Fval <- T2 * (n - p) / (p * (n - 1))
pval <- pf(Fval, p, n - p, lower.tail = FALSE)

cat("Test af H0: mu = [500, 50, 30]\n",
    "Stikprøvestørrelse (n):", n, "\n",
    "Antal variable (p):", p, "\n",
    "Stikprøvegennemsnit:", x.mean, "\n",
    "T^2 statistik:", T2, "\n",
    "F-værdi:", Fval, "\n",
    "p-værdi:", pval, "\n",
    "Signifikansniveau (alpha):", 0.05, "\n",
    "Beslutning:", ifelse(pval < 0.05, "Forkast H0", "Forkast ikke H0"), "\n")
## Test af H0: mu = [500, 50, 30]
##  Stikprøvestørrelse (n): 87 
##  Antal variable (p): 3 
##  Stikprøvegennemsnit: 526.5862 54.68966 25.12644 
##  T^2 statistik: 223.3102 
##  F-værdi: 72.70564 
##  p-værdi: 0.00000000000000000000002828097 
##  Signifikansniveau (alpha): 0.05 
##  Beslutning: Forkast H0

Baseret på testen kan vi konkludere, at der er tilstrækkelig evidens til at forkaste nulhypotesen \(H_0: \boldsymbol{\mu}' = [500, 50, 30]\) på \(\alpha = .05\) signifikansniveau. P-værdien er mindre end signifikansniveauet, hvilket betyder, at forskellene mellem de observerede stikprøvegennemsnit og de forventede middelværdier er statistisk signifikante.

Dette indikerer, at der er grund til at tro, at gruppen af studerende repræsenteret ved scorerne i Tabel 5.2 scorer anderledes end de tusinder af college-studerende over de sidste 10 år, når vi ser på alle tre variable samtidigt. Selvom stikprøvegennemsnittene (526.5, 54.59, 25.13) er “relativt” tæt på de forventede værdier (500, 50, 30), kan forskellen ikke forklares ved tilfældig variation.

Bonusopgave 5.18.b

Bestem længderne og retningerne for akserne af 95% konfidensellipsoide for \(\boldsymbol{\mu}\).

# Beregner den kritiske værdi for 95% konfidensellipsoide
T2cv <- (n - 1) * p / (n - p) * qf(0.95, p, n - p)

# Beregner egenværdier og egenvektorer af kovariansmatricen
r <- eigen(S)

# Beregner længderne af akserne
axislengths <- sqrt(r$values * T2cv / n)

cat("Kritisk værdi af T^2:", T2cv, "\n")
## Kritisk værdi af T^2: 8.333483
# Akselængder
print(axislengths)
## [1] 23.729998  2.472768  1.182500
# Retninger af akserne (egenvektorer)
print(r$vectors)
##            [,1]         [,2]         [,3]
## [1,] 0.99390539  0.103731534 -0.037307396
## [2,] 0.10344339 -0.994589227 -0.009577815
## [3,] 0.03809906 -0.005660238  0.999257936

Konfidensellipsoiden for middelværdivektoren \(\boldsymbol{\mu}\) er karakteriseret ved sine akser. Længderne af disse akser er beregnet som \(\sqrt{\lambda_i \cdot T^2_{cv} / n}\), hvor \(\lambda_i\) er egenværdierne af kovariansmatricen, og \(T^2_{cv}\) er den kritiske værdi for T²-fordelingen.

Bonusopgave 5.18.c

Konstruer Q-Q plots fra marginalfordelingerne af social science and history, verbal, og science scores. Konstruer også de tre mulige spredningsdiagrammer fra parrene af observationer på forskellige variabler. Ser disse data ud til at være normalfordelte? Diskuter.

# Q-Q plots for de marginale fordelinger
par(mfrow=c(1,3))
for (i in 1:p) {
  qqnorm(tbl[,i], main = paste("Q-Q Plot for", colnames(tbl)[i]))
  qqline(tbl[,i])
}

# Spredningsdiagrammer for par af variable
par(mfrow=c(1,3))
for (i in 1:(p-1)) { 
  for (j in (i+1):p) {
    plot(tbl[,i], tbl[,j], 
         main = paste(colnames(tbl)[i], "vs", colnames(tbl)[j]),
         xlab = colnames(tbl)[i], ylab = colnames(tbl)[j])
  }
}

# chi-i-anden plot for multivariat normalitet
par(mfrow=c(1,1))
Sx <- cov(tbl)
D2 <- mahalanobis(tbl, colMeans(tbl), Sx)
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis-afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

Baseret på de grafiske analyser kan vi vurdere, om dataene følger en multivariat normalfordeling:

  1. Q-Q plots for de marginale fordelinger:
    • Social Science and History: Punkterne følger den rette linje rimeligt godt, med små afvigelser i halen.
    • Verbal: Punkterne følger den rette linje meget tæt, men med en “systematisk” afvigelse i halen.
    • Science: Generelt en fin tilpasning.
  2. Spredningsdiagrammer:
    • Punkterne i spredningsdiagrammerne viser ingen tydelige mønstre, der ville indikere ikke-normalitet.
    • Der ser ud til at være positive korrelationer mellem variablerne, hvilket er forventeligt for test-scores.
  3. chi-i-anden plot:
    • Punkterne følger den rette linje godt for mindre værdier, men med tydelige afvigelser i den højre hale.

Konklusion: Overordnet set viser de grafiske analyser, at dataene ikke umiddelbart følger en multivariat normalfordeling. Der er afvigelser, især i halerne af marginalfordelingen for Verbal og i chi-i-anden plottet. Især afvigelserne i chi-i-anden plottet er alvorligt nok til at afvise antagelsen om multivariat normalitet. Derfor er brugen af Hotelling’s T² test og konstruktionen af konfidensellipsoiden ikke nødvendigvis berettiget for disse data.

Bonusopgave 5.19

Målinger af \(X_1 =\) stivhed og \(X_2 =\) bøjningsstyrke for en stikprøve på \(n = 30\) stykker af en bestemt kvalitet af tømmer er givet i Tabel 5.11. Enhederne er pounds/(inches)².

Bonusopgave 5.19.a

Konstruer og skitser en 95% konfidensellipse for parret \([\mu_1, \mu_2]'\), hvor \(\mu_1 = E(X_1)\) og \(\mu_2 = E(X_2)\).

# Indlæser data
data("T5-11")

# Beregn nødvendige statistikker
n <- nrow(tbl)
p <- ncol(tbl)
x.mean <- colMeans(tbl)
S <- cov(tbl)

# Definer et grid for at plotte konfidensellipsen
x <- seq(x.mean[1] - sqrt(S[1, 1]), x.mean[1] + sqrt(S[1, 1]), length.out = 100)
y <- seq(x.mean[2] - sqrt(S[2, 2]), x.mean[2] + sqrt(S[2, 2]), length.out = 100)
nx <- length(x)
f <- matrix(0, nx, nx)

# Beregn Mahalanobis-afstande for hvert punkt i grid'et
for (i in 1:nx) {
  for (j in 1:nx) {
    mu <- c(x[i], y[j])
    f[i, j] <- n * t(x.mean - mu) %*% solve(S) %*% (x.mean - mu)
  }
}

# Beregn den kritiske værdi for 95% konfidensellipse
T2cv <- (n - 1) * p / (n - p) * qf(0.95, p, n - p)

# Plot konfidensellipsen
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 2, 
        xlim = c(min(tbl[,1]), max(tbl[,1])), ylim = c(min(tbl[,2]), max(tbl[,2])),
        xlab = "Stivhed (pounds/inches²)", ylab = "Bøjningsstyrke (pounds/inches²)",
        main = "95% Konfidensellipse for middelværdier af stivhed og bøjningsstyrke")

# Tilføj datapunkter og middelværdi
points(tbl, pch = 1)
points(x.mean[1], x.mean[2], pch = 16, col = "red")

# Vis statistikker
cat("95% konfidensellipse for middelværdivektoren [μ₁, μ₂]':\n",
    "Stikprøvestørrelse (n):", n, "\n",
    "Antal variable (p):", p, "\n",
    "Stikprøvegennemsnit af stivhed (x̄₁):", x.mean[1], "\n",
    "Stikprøvegennemsnit af bøjningsstyrke (x̄₂):", x.mean[2], "\n",
    "Kritisk værdi af T²:", T2cv, "\n")
## 95% konfidensellipse for middelværdivektoren [μ₁, μ₂]':
##  Stikprøvestørrelse (n): 30 
##  Antal variable (p): 2 
##  Stikprøvegennemsnit af stivhed (x̄₁): 1860.5 
##  Stikprøvegennemsnit af bøjningsstyrke (x̄₂): 8354.133 
##  Kritisk værdi af T²: 6.91937
# Beregn egenvektorer og egenværdier for at vise akserne
eigSol <- eigen(S)
axislengths <- sqrt(eigSol$values * T2cv / n)

# Akselængder
cat("Akse 1:", axislengths[1], "\n",
    "Akse 2:", axislengths[2], "\n")
## Akse 1: 901.6522 
##  Akse 2: 140.5119

Figuren viser en 95% konfidensellipse for middelværdivektoren \([\mu_1, \mu_2]'\), hvor \(\mu_1\) er middelværdien af stivhed og \(\mu_2\) er middelværdien af bøjningsstyrke. Ellipsen er konstrueret baseret på 30 observationer (vist som sorte cirkler), og det røde punkt repræsenterer stikprøvegennemsnittet.

Ellipsens form og orientering afspejler kovariansstrukturen mellem de to variable. Den er mere udstrakt langs akse 1, hvilket indikerer større usikkerhed i denne retning. Akserne af ellipsen svarer til egenvektorerne af kovariansmatricen, og deres længder er relateret til egenværdierne.

Bonusopgave 5.19.b

Antag at \(\mu_{10} = 2000\) og \(\mu_{20} = 10,000\) repræsenterer “typiske” værdier for stivhed og bøjningsstyrke. Givet resultatet i (a), er dataene i Tabel 5.11 konsistente med disse værdier? Forklar.

# Definer de "typiske" værdier
mu0 <- c(2000, 10000)

# Plot konfidensellipsen
contour(x, y, f, levels = T2cv, drawlabels = FALSE, lwd = 2, 
        xlim = c(min(tbl[,1]), max(tbl[,1])), ylim = c(min(tbl[,2]), max(tbl[,2])),
        xlab = "Stivhed (pounds/inches²)", ylab = "Bøjningsstyrke (pounds/inches²)",
        main = "95% Konfidensellipse for middelværdier af stivhed og bøjningsstyrke")

# Tilføj datapunkter og middelværdi
points(tbl, pch = 1)
points(x.mean[1], x.mean[2], pch = 16, col = "red")

# Marker punktet i plottet
points(mu0[1], mu0[2], pch = "X", col = "blue", cex = 1.5)

# Beregn den kvadrerede Mahalanobis-afstand for de typiske værdier
D2 <- n * t(x.mean - mu0) %*% solve(S) %*% (x.mean - mu0)

# Sammenlign med den kritiske værdi
cat("Test af om de typiske værdier er konsistente med data:\n",
    "Typiske værdier: [μ₁₀, μ₂₀]' = [", mu0[1], ",", mu0[2], "]'\n",
    "Kvadreret Mahalanobis-afstand:", D2, "\n",
    "Kritisk værdi (T²):", T2cv, "\n",
    "Er værdierne inden for konfidensellipsen?", D2 < T2cv, "\n")
## Test af om de typiske værdier er konsistente med data:
##  Typiske værdier: [μ₁₀, μ₂₀]' = [ 2000 , 10000 ]'
##  Kvadreret Mahalanobis-afstand: 23.64778 
##  Kritisk værdi (T²): 6.91937 
##  Er værdierne inden for konfidensellipsen? FALSE

De “typiske” værdier \(\mu_{10} = 2000\) og \(\mu_{20} = 10,000\) er markeret på plottet med et blåt X. For at afgøre, om disse værdier er konsistente med vores data, sammenligner vi den kvadrerede Mahalanobis-afstand for disse værdier med den kritiske værdi for T²-statistikken.

Da den kvadrerede Mahalanobis-afstand er større end den kritiske værdi, ligger punktet (2000, 10000) udenfor 95% konfidensellipsen. Dette betyder, at de “typiske” værdier IKKE er konsistente med dataene i tabellen, og vi kan forkaste hypotesen om, at den sande middelværdivektor er \([\mu_1, \mu_2]' = [2000, 10000]'\).

Med andre ord, baseret på vores data er det plausibelt, at den gennemsnitlige stivhed er 2000 pounds/inches² og den gennemsnitlige bøjningsstyrke er 10,000 pounds/inches² for denne kvalitet af tømmer.

Bonusopgave 5.19.c

Er den bivariate normalfordeling en brugbar populationsmodel? Forklar med henvisning til Q-Q plots og et spredningsdiagram.

# Q-Q plots for de marginale fordelinger
par(mfrow=c(1,2))
for (i in 1:p) {
  qqnorm(tbl[,i], main = paste("Q-Q Plot for", colnames(tbl)[i]))
  qqline(tbl[,i])
}

# Spredningsdiagram
par(mfrow=c(1,1))
plot(tbl, main = "Spredningsdiagram af stivhed mod bøjningsstyrke", 
     xlab = "Stivhed (pounds/inches²)", ylab = "Bøjningsstyrke (pounds/inches²)")

# chi-i-anden plot for multivariat normalitet
Sx <- cov(tbl)
D2 <- mahalanobis(tbl, colMeans(tbl), Sx)
qqplot(qchisq(ppoints(n, a = 0.5), df = p), D2,
       ylab = "Mahalanobis-afstande",
       xlab = bquote("Kvantiler af " ~ chi[.(p)]^2),
       main = bquote("Q-Q plot af Mahalanobis" * ~ D^2 * 
                       " vs. kvantiler af" * ~ chi[.(p)]^2))
abline(0, 1)

Baseret på de grafiske analyser kan vi vurdere, om den bivariate normalfordeling er en brugbar model for disse data:

Konklusion: Overordnet set viser de grafiske analyser, at den bivariate normalfordeling er en brugbar populationsmodel for disse tømmerdata. De marginale fordelinger er tilnærmelsesvis normale, spredningsdiagrammet viser et elliptisk mønster, og chi-i-anden plottet bekræfter, at de multivariate afstande følger den forventede fordeling. Der er små afvigelser i halerne, men disse er ikke alvorlige nok til at afvise antagelsen om bivariat normalitet. Derfor er brugen af T²-baserede konfidensregioner og -tests berettiget for disse data.

Bonusopgave 5.30

Referer til dataene om energiforbrug i Opgave 3.18.

Bonusopgave 5.30.a

Bestem large-sample 95% Bonferroni-konfidensintervaller for det gennemsnitlige forbrug af hver af de fire typer, summen af de fire, og forskellen, petroleum minus naturgas.

# Kovariansmatricen er givet som en nedre triangulær matrix
S <- diag(4)
S[lower.tri(S, diag = TRUE)] <- c(0.856, 0.635, 0.173, 0.096, 0.568, 0.128, 0.067, 0.171, 0.039, 0.043) 
S <- S + t(S) - diag(diag(S))
print(S)
##       [,1]  [,2]  [,3]  [,4]
## [1,] 0.856 0.635 0.173 0.096
## [2,] 0.635 0.568 0.128 0.067
## [3,] 0.173 0.128 0.171 0.039
## [4,] 0.096 0.067 0.039 0.043
# Middelværdivektoren
x.mean <- c(0.766, 0.508, 0.438, 0.161)

# Stikprøvestørrelse og antal variable
n <- 50
p <- length(x.mean)

# Antal sammenligninger (4 individuelle gennemsnit + summen + forskellen)
m <- 6

# Beregner den kritiske værdi for Bonferroni-justeringen
ncv <- qnorm(1 - 0.05 / (2 * m))

# Beregner Bonferroni-konfidensintervaller for de individuelle gennemsnit
nBimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- ncv * sqrt(S[i, i] / n)
  nBimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(nBimu) <- c("Kul", "Naturgas", "Petroleum", "Elektricitet")
colnames(nBimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregner Bonferroni-konfidensinterval for summen
csum <- c(1, 1, 1, 1)
xmsum <- as.vector(t(csum) %*% x.mean)
Ssum <- as.vector(t(csum) %*% S %*% csum)
cisum <- c(xmsum, xmsum - ncv * sqrt(Ssum / n), xmsum + ncv * sqrt(Ssum / n))
names(cisum) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregner Bonferroni-konfidensinterval for forskellen (petroleum minus naturgas)
cdif <- c(0, -1, 1, 0)
xmdif <- as.vector(t(cdif) %*% x.mean)
Sdif <- as.vector(t(cdif) %*% S %*% cdif)
cidif <- c(xmdif, xmdif - ncv * sqrt(Sdif / n), xmdif + ncv * sqrt(Sdif / n))
names(cidif) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# 95% Bonferroni-konfidensintervaller for gennemsnitsforbrug af de fire energityper
print(nBimu)
##              Middelværdi Nedre grænse Øvre grænse
## Kul                0.766   0.42080140   1.1111986
## Naturgas           0.508   0.22680583   0.7891942
## Petroleum          0.438   0.28371269   0.5922873
## Elektricitet       0.161   0.08363111   0.2383689
# 95% Bonferroni-konfidensinterval for summen af de fire energityper
print(cisum)
##  Middelværdi Nedre grænse  Øvre grænse 
##     1.873000     1.134854     2.611146
# 95% Bonferroni-konfidensinterval for forskellen (petroleum minus naturgas)
print(cidif)
##  Middelværdi Nedre grænse  Øvre grænse 
##   -0.0700000   -0.3293019    0.1893019

Bonusopgave 5.30.b

Bestem large-sample 95% simultane \(T^2\)-intervaller for det gennemsnitlige forbrug af hver af de fire typer, summen af de fire, og forskellen, petroleum minus naturgas. Sammenlign med dine resultater fra Del a.

# Beregner den kritiske værdi for T^2
Xcv <- qchisq(1 - 0.05, p)

# Beregner T^2 konfidensintervaller for de individuelle gennemsnit
nTimu <- matrix(0, p, 3)
for (i in 1:p) {
  err <- sqrt(Xcv) * sqrt(S[i, i] / n)
  nTimu[i,] <- c(x.mean[i], x.mean[i] - err, x.mean[i] + err)
}
rownames(nTimu) <- c("Kul", "Naturgas", "Petroleum", "Elektricitet")
colnames(nTimu) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregner T^2 konfidensinterval for summen
ciTsum <- c(xmsum, xmsum - sqrt(Xcv) * sqrt(Ssum / n), xmsum + sqrt(Xcv) * sqrt(Ssum / n))
names(ciTsum) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Beregner T^2 konfidensinterval for forskellen
ciTdif <- c(xmdif, xmdif - sqrt(Xcv) * sqrt(Sdif / n), xmdif + sqrt(Xcv) * sqrt(Sdif / n))
names(ciTdif) <- c("Middelværdi", "Nedre grænse", "Øvre grænse")

# Viser resultaterne
# 95% T^2 konfidensintervaller for gennemsnitsforbrug af de fire energityper
print(nTimu)
##              Middelværdi Nedre grænse Øvre grænse
## Kul                0.766   0.36297404   1.1690260
## Naturgas           0.508   0.17970044   0.8362996
## Petroleum          0.438   0.25786662   0.6181334
## Elektricitet       0.161   0.07067034   0.2513297
# 95% T^2 konfidensinterval for summen af de fire energityper
print(ciTsum)
##  Middelværdi Nedre grænse  Øvre grænse 
##       1.8730       1.0112       2.7348
# 95% T^2 konfidensinterval for forskellen (petroleum minus naturgas)
print(ciTdif)
##  Middelværdi Nedre grænse  Øvre grænse 
##   -0.0700000   -0.3727399    0.2327399